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Abstract. Cells can often choose among several stably heritable phenotypes. 
Examples are the expression of genes in eukaryotic cells where long chromosomal 
q ' regions can adopt persistent and heritable silenced or active states, that may be 

associated with positive feedback in dynamic modification of nucleosomes. We 
generalize this mechanism in terms of bistability associated with valleys in an epigenetic 
landscape. A transfer matrix method was used to rigorously follow the system through 
the disruptive process of cell division. This combined treatment of noisy dynamics both 
between and during cell division provides an efficient way to calculate the stability of 
■ alternative states in a broad range of epigenetic systems. 
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Introduction 

X 

Cells carry information handed down from their ancestors and are able to pass on 
information to their descendants. In many cases this "memory" is epigenetic, that 
is, not stored in the DNA sequence, allowing cells with identical DNA to maintain 
distinct functional identities. Epigenetic cell memory implies alternative states of gene 
expression that are stable over time and are inherited through cell division. 

A proposed mechanism for epigenetic cell memory invokes positive feedback loops 
in nucleosome modification [H El El HI EJ E] • Positive feedback is a mechanism seen in 
many other regulatory systems where for example production of a regulatory protein 
activates its own production, or more robustly where two mutual repressors act strongly 
enough to prevent co-expression. A complementary view on cell memory is that of an 
epigenetic landscape [3 [8] , where the state of a cell develops on some potential energy 
surface, and a state is maintained when the cell is caught at a particular valley for a 
long time. 
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In this paper we develop an epigenetic landscape formalism for cell memory by 
positive feedback in nucleosome modification. Instead of viewing cell differentiation as 
the 'pushing' of a cell over a fixed landscape [TJ [8] , our approach suggests that cell fate 
could be controlled by changing the landscape. 

Models 

Inspired by the mating type switch in S. pombe [9] we introduced a model for bistability 
by positive feedback in nucleosome modification [I] . The model had one parameter, the 
positive feedback to noise ratio F, and modeled the dynamics of a system consisting 
of N nucleosomes where each could be in one of 3 states, modified, unmodified and 
anti-modified. 

Here we introduce a simpler version of this modification system, in which there are 
only two chemical states of each nucleosome. We term these modified (M) and anti- 
modified (A) to indicate their mutual exclusivity, with the A nucleosome carrying either 
a different chemical modification or no modification (Fig. CD). Each nucleosome type 
recruits a modifying enzyme that converts the other type to its own type. Our results 
from [1] demonstrated that robust bistability requires an effective cooperativity in the 
recruitment process. Cooperativity could here be included directly by requiring that 2 
local nucleosomes with the same modification, e.g. M, are needed to make an A— >M 
conversion, as described in Fig. [U Note that this model is not only just a simplified 
version of 3-state model but also has parallels with the mating-type silencing system in 
S. cerevisiae, where one typically considers acetylated and non-acetylated nucleosomes 

piainra. 

Suppose M nucleosomes are in the M-state at time t. We now express the 
development of the fraction m = M/N sites in the M-state. Denoting dm = 1/N, 
we have 

dm 

—jj- = (R + (m) — i?_(m)) • dm + noise (1) 

with the rate that the system with fraction m of M-sites gets one more (or one less) 
M-sites being R + (m) (or i?_(m)) and with noise having zero mean and being associated 
with the randomness of processes in a finite system. The rates are given by 

R + (m) = a(l — m)m 2 + (1 — a)(l — m), 

R_{m) = am{l — m) 2 + (1 — a)m. (2) 

Here, the first term is the nucleosome recruitment; in the case of R + (m), the recruitment 
occurs with probability a, and it must involve 2 M-sites (probability proportional to m 2 ) 
and must change the modification on an A-site (probability proportional to (1 — m)). 
The second term is the noise effect which is proportional to (1 — a), where a nucleosome 
can become M by random conversion of from an A state (probability proportional to 
(1 — m)), and vice versa. This noise represents all events from the cell that are not 
associated to the direct recruitment processes from other nucleosomes within our N 
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nucleosome system. One should be aware that there is another level of noise in our 
stochastic description, represented by the noise in eq.flT]) which is the noise associated 
with the stochasticity of the molecular processes. 

The ratio of recruitment (or the positive feedback) to the noise, F = a/(l — a) is 
the parameter of the model. The Langevin equation for m is then 

—jj- = — (2m — l)[m(l — m) — 1/F] + noise, (3) 

implicitly implying that there is bistability when F > 4 provided that the noise term 
is small (system size large). 

To analyze the time development of the distribution of m more carefully we 
reformulate the model in terms of a master equation for the probability P(m, t) as 

d 

—P(m, t) = R_(m + dm)P{m + dm, t) + R + {m — dm)P{m — dm, t) 
at 

-[R + (m)+R-(m)]P(m,t). (4) 

Results 

Epigenetic landscape generated from a positive feedback system 

We extract "the potential landscape" in m-space by comparing eq. (|4j) with the generic 
1-dimensional Fokker-Planck equation for diffusion of a particle in a potential U(m): 



d _ dJ d 



dt dm dm 

_d_ 
dm 



-fi(m) dU ^ P - d ( D ( m ) P ) 
dm dm 

. dVjjn) dP 

^ m) ^m- P D{m) dm- 



(5) 



Here, fi(m) is the mobility and D(m) quantifies the stochastic motion in terms of an m- 
dependent diffusion coefficient. In the last step, V(m) represents an effective potential 
that includes both drift and noise events, defined as ^- = 4^- + H dX °^ D \ 

dm am dm 

Expanding eq. 02]) with eq. (jl]) to second order in dm = 1/N and comparing it with 
eq. ([3]), we find the drift (^f ) = p,(m)(dU / dm) , the effective potential V(m), diffusion 
D(m) and mobility fi(m) as the following: 

dm \ rv 

v =i± (2m-l)[m(l-m)-l/F] (6) 

= 2JVm(l - m) + M - — J ln[Fm(l - m) + 1] (7) 

. , . , m(l — m) + l/F , , 

D(m) =fM(m) = a 1 ^ 1 . (8) 

Here, the first equation could have been obtained directly from the Langevin equation 
(J3]). From these expressions we again can see that there is the critical recruitment to 
noise ratio F = a/(l — a), with m — 1/2 being an unstable fixed point for F > 4. 

Figure |2] shows (dm/dt), V(m) and the steady state distribution Po(m) for F = 3 
and F = 12, thereby illustrating mono-stable and bistable systems. Also note that the 
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analytic results fit the stochastic simulation, with a deviation that scales as 1/N with 
increased system size (not shown). Figure [3] shows how the epigenetic landscape changes 
gradually as F increases: from a single steep valley, through an almost equipotential 
'river plain', to two valleys. Because F depends on protein concentrations and affinities, 
the shape of the epigenetic landscape is under biological control. 

One can repeat these calculations for a model where recruitment is not requiring 
the cooperative action of two nucleosomes (The second order terms in eq. ([2]) should 
then be replaced by first order). In that case one never obtains more than one stable 
fixed point, confirming that bistability indeed requires cooperativity [UEj. 

The potential V(m) effectively describes the effective force on m from the combined 
effect of recruitment and noise events. Thus a large positive gradient in V(m) means that 
nucleosomes in average will tend to lose their m modification. A potential minimum, on 
the other hand, means that recruitment processes and noise events balance such that 
the number of modified states typically stays around this minimum. In this way our 
potential V(m) plays the role of a epigenetic landscape in the Waddington sense [7j. 
In particular, the valleys and hills of this landscape can be viewed as the metastable 
epigenetic states and the barriers between them. We will use this analogy to calculate 
first the probability for stochastic switching between such states, and subsequently we 
will discuss how one may alter the landscape by modifying the recruitment processes 
that define the landscapes. A modification that was also envisioned by strings in the 
Waddington landscape!?]. 



Stability of a macroscopic state 

Now we quantify the stability of a macroscopic state by the average number of attempted 
updates per nucleosome before the full system switches for the first time to the alternate 
epigenetic state. Using D = fi from eq. (ED, we in analogy to Kramers [13] rewrite 



J 



-D 



dm dm 



d 

-Dexp(-V)— LPexp(V)] 
dm 



(9) 



and use the quasi- stationary approximation (i.e., the current J is constant) to write the 
flux for going from an A-state (the potential minimum at m — mA ~ 0) to an M-state 
(the potential minimum at m = mu ~ 1): 



J 



[Pexp{V)]? M 



C™(l/D(m))exp(V(m))dm- 



(10) 



Using a Gaussian approximation (i.e. V(m) harmonic around both the initial state A 
and the transition state T with m = mx = 1/2 and the initial distribution for P(m,t) 
around the state A), we obtain the average life time of an epigenetic state r as 



-L^ATrN^exp[V T -V A ] 



(11) 
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for large N and F, where Vr = Vim = 1/2), and Va = V(rriA) is the potential 
minimum for the A state (The detailed calculation is given in the Appendix). Figure H] 
demonstrates that eq. (TTT]) reproduces stochastic simulations. However, when pushing 
towards very small N, there is a tendency that the continuous description deviate from 
the stochastic result. Thus, for iV of order 10 or below, we recommend a stochastic 
simulation. 

Equation ( II ip can also be used to obtain an interesting prediction from our 
model. Using the expression for the potential V from eq. (J7J) for large N, we see that 
V(m — ~) ~ Nf(F) with a function f(F) independent of N, and thus that stability 
scales exponentially with N, i.e., r oc 

Effect of cell divisions 

Epigenetic states are capable of being inherited across cell divisions. This can give 
difficulties for stability of the states |4], particularly for 2-state systems [14] At cell 
division the genome is duplicated, and following refs. [U [T5] we assume that the 
resident nucleosomes are partitioned randomly between the daughter strands. The 
vacant positions are filled by new randomly selected nucleosomes where half are in 
M and half in A-state. We accordingly supplement our model above with cell divisions 
at certain fixed time intervals. This cell generation time is measured in units of the 
number of attempted nucleosome updates per nucleosome. 

Whereas the potential landscape between cell divisions drives the system toward 
one of the epigenetic states, the randomization at cell divisions brings the system closer 
to the top of the potential barrier in the epigenetic landscape. 

Consider that before cell division the system is in a state with M& = nib x N 
nucleosomes in the M-state and the remaining nucleosomes in the A state. Cell division 
results in the distribution of number of M-state nucleosomes M a 



where the sum runs over all the ways of getting from M& to M a by selecting M < 
min(M a , Mb) nucleosomes in the M-state and A < min(iV — M a , N — Mb) nucleosomes 
in the A state to be transferred directly at the cell division. 

Between cell divisions, the system evolves by a stochastic sequence of single 
nucleosome exchanges that can be described by motion in the epigenetic landscape. 
The stochastic change M — » M ± 1 can be followed by the master equation above, and 
the system evolves towards a steady state distribution with two well separated peaks, 
see Fig. [2j 

To combine the gradual development in a well defined epigenetic landscape between 
cell division, with the sudden reshuffling at cell divisions we express the gradual 
development in terms of matrix operations. The matrix G that correspond to eq. (J2]) 
is only non-zero at the diagonal and off-diagonal elements: G(M,M + 1) = R + (m), 



G(M,M - 1) = R-{m) and G(M,M) = 1 - G(M,M - 1) - G(M,M + 1) where 



D{M a ,M b ) = Y,H 



( M bM) ( N — M b A) (N - M- AM a - M) 



(12) 



M A 



2{2N-M-A) 
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m = M/N and R±(m) is from eq. (j2J). The probability distribution evolves according to 
P(M) -> P(M') = £ M G(M, M')P(M) for each update of a nucleosome in the system. 

In Fig. [5] we show the time evolution of the probability distribution from one cell 
division to the next for a system with iV = 60, F = 10 and m ~ 0. In panel A, 
the simulation uses a generation length of 20 updates per nucleosome. The evolution 
starts just after a cell division, where a randomization (using eq. [12]) is followed by 
a drift imposed by the epigenetic landscape. Just before the next cell division, one 
resets P(M) = for M > N/2, and renormalizes the distribution. We see that at cell 
division a small fraction of cells reach large m values, and over the next ~ 10 updates 
per nucleosome can move to m ~ 1. After around 10 updates, the P(m) distribution 
reaches a quasi-steady state, reflecting that from then on a very small flux goes over 
the barrier. Thus after 10 updates the likelihood of further transitions between the two 
epigenetic states can be ignored (for F > 10). 

The time evolution in Fig. [5] also illustrates that transitions are entirely dominated 
by the noise at cell divisions, at least for large enough F. In fact by stochastic simulation 
we have verified that switches only occur when the stochastic partitioning in a division 
brings the system close to the transition state, m ~ 1/2. 

An entire cell generation with cell division is described by the matrix 

C(M, M") = ]T G& N \M, M')D(M', M") (13) 

M' 

where g is the number of single nucleosome updates per cell generation. Iterating the 
updating process with renormalization of the distribution, we obtain the probability 
distribution P(M') shortly before a cell division, and estimate the average number of 
generations needed before escape (escape time) as n e = (J2m'>n/2 P{M')\ in the panel 
B. This escape time is the average time it takes for a system to switch from one state 
to another. From Figure [5] we see that when g > 10, the escape time does not depend 
on the value of g for F > 10 because of the small escaping rate after 10 updates. 

Epigenetic landscape with regulated tilt 

Finally, we consider the case where the modification and anti-modification are not 
symmetric but one of the effects is stronger than the other. This asymmetry of 
modification could be under biological control in a real system by changing the 
concentration of modifying enzymes or by recruiting such enzymes by transcription 
factors [Bj. 

One of the simplest ways to include such asymmetry is to modify the transition 
rate eqs. (|2D into 

R + (m) = 2(1 — r))a(l — m)m 2 + (1 — ot)(l — to), 

i?_(m) = 2r]am(l — m) 2 + (1 — a)m. (14) 

Here, a new parameter 77, defined in the range < 77 < 1, sets the relative strength of 
modification versus anti-modification, and gives the symmetric case when 77 = 1/2. 
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Figure shows how the potential landscape can be manipulated by changing 
the relative strengths of the two recruitment processes, rj. This potential is directly 
defined as V(m) = — In P Q (m), where Po(m) is defined as the steady state probability 
distribution of m. For simplicity the cell division is not taken into account here. For 
small r], the valley at the M-state is much deeper than that of the A-state, and the 
landscape dramatically changes as rj increases to give a landscape where the A-state has 
a deeper valley than the M-state. 

Discussion 

Positive feedback in nucleosome modification is a powerful mechanism to maintain 
a dynamic bistable system, even with destabilizing factors such as cell-division. 
Here we demonstrated how positive feedback in itself can be reformulated into an 
epigenetic landscape with peaks and valleys that reflects the underlying balance between 
feedback and noise. As long as movements are small, dominated by single nucleosome 
modifications, the movement in the landscape can be fully modeled by a Langevin or 
Fokker-Planck equation with a first escape time calculated in analogy to Kramers. When 
stochastic events are large, as during cell divisions, a transfer matrix method allowed us 
to extend the Fokker-Planck formalism and thereby to set a minimum timescale for the 
dynamics of a robust positive feedback. Finally, we studied how the landscape could be 
'tilted' by asymmetry in the nucleosome modification reactions. 

We expect that the transfer matrix method can be extended relatively easily to 
include other nucleosome modification schemes, for example the 3-state model of p[]. 
This approach is more powerful than the mean field approach [13] in the sense that it 
allows one to explore the probability distribution. The method can also deal with cases 
where the epigenetic system does not conform to a potential energy surface. Noise 
associated with cell divisions can for example be included in systems with double 
negative feedback between repressors, such as the Cl-Cro feedback loop in the lysis- 
lysogeny switch of phage lambda [16J. 

Our results are consistent with recent observations in mammalian cells in which 
increased cell division rates accelerated stochastic transitions between epigenetic states 

Epigenetic landscapes present a particularly appealing way to discuss multi-stability 
of expression states in living systems. The presented coupling between positive feedback 
and the possibility for a drift in a landscape, may be useful for understanding cases 
where bistable decisions are delayed, as often seems to be the case in development. 
Some epigenetic landscapes may define the activity of transcription factors that act as 
histone modifying complexes, and thereby subsequently define the input parameters for 
other landscapes further down along a developmental pathway. Thereby understanding 
of epigenetic stability and regulated tilting of landscapes may speak to large classes of 
coupled switch systems. 
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Conclusion 

This paper explored theoretical implications of epigenetics as a dynamic phenomenon, 
where alternate states of gene expression are selected and maintained over multiple 
generations by on going dynamic processes. Our approach builds on the assumption 
that the epigenetic states are maintained through a positive feedback where nucleosomes 
of a certain kind recruit enzymes which in turn convert other nucleosomes to the same 
kind. By introducing a minimal model for such a dynamic system we demonstrated 
that the on-going nucleosome updating rates only need to be g ~ 10 updates per 
nucleosome per generation to provide robust maintenance. This result is closely linked 
to the convergence of the P(M) distribution after cell division (Fig. 4), which we find 
to be rather insensitive to the value of F. Thus we expect that any model working with 
positive feedback as a key maintenance factor, in principle would work with a moderate 
number of updates per generation. 

It is important to notice that the stability of the epigenetic state depends on the 
update rules at the cell-divistion especially in the two-state model. In our model, on 
average half of the nucleosomes come from parental DNA, and the rest of them are 
either in the modified or the unmodified state with equal probability, thus effectively 
75% nucleosome are in the same epigenetic state as the parental DNA. Thus the most 
of the cells are still in the same valley in the epigenetic potential (Fig. 2 g) and come 
back to the original epigenetic state. However, if half of the nucleosomes are replaced 
with unmodified nucleosomes at the cell division in the 2-state model, the system prefers 
unmodified states and one needs extra mechanisms to keep modified state stable [14J. 
It is not clear which mechanism operates at DNA replication in different systems, but 
our formulation is applicable in both cases. More experimental information about the 
nucleosomes inserted after DNA replication will be critical in understanding the stability 
of epigenetic states. 

Our analysis also showed formally that bistability requires cooperativity of 
recruitment, in the sense that eq. ([2]) require more than first order terms in 
order to provide separation between two states. In the model in Dodd [I] the 
recruitment was implemented by requiring that two independent recruitment processes 
were of the same type, whereas we here simply assumed that the two simultaneously 
recruiting nucleosomes are of the same type. Our analogy between system size and 
effective randomness allowed us to show formally that stability of inherited states will 
grow exponentially with the number of nucleosomes in the considered region of the 
chromosome. 

Our analysis opens for understanding development in epigenetic landscapes, in 
terms of positive feedback mechanisms that are linked to each other through expression 
of nucleosome modifying enzymes. 
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Appendix 

Analytical calculation of the stability of a macroscopic state 

Assuming V(m) is harmonic around m = m A and hence the initial distribution P(m 
is a Gaussian around m = m A (and P(rriM) ~ 0), we get 

1 



[P{m) exp(V(m))] 



m M 
m A 



-P(m A )exp(V(m A )) 



2na 2 A 



m=m,A 



From the condition dV(m)/dm 



with l/a 2 A = (d 2 V(m)/dm 2 

©, we have m A = 1/2 -yJ\/A- 1/F+ 1/(2JV), and 

1 AN -F 



m=m,A 



exp(V(m A )),(15) 
and eq. 



o-- 



'd 2 V(m) 
dm 2 



1 



4iV(FiV + 2F- AN) N(F -A) NF* 



(16) 



m=mjs i 

where we assume iV >> 1 and F » 1 (i.e. a ~ 1). 

Approximating V(jn) as harmonic around the transition state T with m = rax 
1/2 and noting V(m) takes the maximum at m — mx, we have 
r m A/ 1 1 



m A 



D(m) 



exp(V(m))dm 



with 



Noting 



we get 



On 



'd 2 V{m) 
dm 2 



D(m T ) 



27T(Ty 



D(m T ) 
1 + F/A 



exp(V^(mT)) 



exp 



(m — tutY 
2o\ 



exp(V(m T )) 



dm 



(17) 



1 + F/A 



m=mj' 



(4 - F)N - 2F N(F - 4) AN 



D(m T ) 



q(l/4 + l/F) 
2N 2 



8N 2 ' 



]N 2 



a 



2na T e v{mT) ■ V2^a A e~ v{mA) 



AtiNJ i exp [V{m T ) - V{m A )\ 



(19) 

(20) 
(21) 
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Figure legends 



Figure 1. Model. Consider a system with N sites. Each site represents a 
nucleosome that can be modified (M) or anti- modified (A). Transitions between these 
two states are in part random, and in part auto-regulated by recruitment of histone 
modifying enzymes by local nucleosomes: At each update, a nucleosome i is either, 
with probability 1 — a, set to an M or an A state randomly. Or with probability a, 
two other nucleosomes are chosen, and if these are in the same state then the state of 
nucleosome i is set to this state. The model is parametrized by the positive feedback 
to noise ratio F = a/(l — a). 
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Figure 2. System properties. Analytical results (solid lines) for a system of size 
N = 60 showing; (A,D): the drift (dm/dt), (B,E): the effective potential V(m) and 
(C,F): the steady state distribution Po(m) in the two regimes: (A,B,C) - F = 3 where 
there is no bistability and (D,E,F) - F = 12 where there is well defined bistability, 
which can be seen in the effective potential (E) with two valleys separated by a hill. 
For (dm/dt) and Po(m) (A,C,D,F), the numerical results are also shown by symbols. 
Comparing the V(m) (B,E) and Po(m) (C,F) notice that ln(P ( m )) ~ —V(m) as 
expected in the steady state from requiring that J = in eq. (j9]). 
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Figure 3. An epigenetic landscape generated from a positive feedback 
system. Here the effective potential V(m) from eq. ([7]) is plotted as a function of 
F with fixed N = 60. The landscape changes gradually as F increases, from a single 
steep valley, through an almost equipotential 'river plain', to two valleys. This change 
is associated with stronger recruitment processes at larger F values. 
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Figure 4. Stability of a macroscopic state. The first passage time, or the average 
life time of an epigenetic state r, is shown as a function of N (A) and F (B). Symbols 
show numerical results and lines show analytical results from eq. Time is counted 
in units of attempted nucleosome updates per nucleosome. 
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Figure 5. The switching of the macroscopic state via randomization at 
the cell division. (A) The averaged development of the M distribution over one 
cell generation, with "time" measured in attempted updates per nucleosome. The 
simulation uses F = 10 and a generation length of 20 updates per nucleosome. The 
evolution starts just after a cell division, where a randomization (using eq. 1121) is 
followed by a drift imposed by the epigenetic landscape. After about 10 updates one 
sees that the P{M) reaches a nearly stationary distribution, where a fixed fraction has 
switched to the alternate state. Just before the next cell division one resets P(M) = 
for M > N/2, and renormalizes the distribution. Iterating this process, panel B shows 
the average number of generations needed before escape (see escape time discussion in 
the text). A direct Monte-Carlo simulation result with a generation time of 30 is also 
shown by symbols. 
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Figure 6. Epigenetic landscape with asymmetry. The figure shows V(m) = 
— ln(P (w)) obtained from the model with asymmetric transition rates eqs. (I14[) using 
the transfer matrix method. For simplicity the cell division is not taken into account 
here. The figure shows how the landscape gradually changes from monostable, to 
bistable and then again to monostable as one changes the parameter eta. This 
"asymmetry" parameter r\ would in principle be under biological control through 
changes of concentrations of modifying enzymes associated with the recruitment 
processes. The total recruitment to noise ratio is fixed at F = 20. The figure illustrate 
a dramatic tilting of the epigenetic landscape by a moderate change in nucleosome 
modification rates. 



